\section{Графы де Брюина и постановка задачи}
\subsection{Основные понятия}

\textit{Риды} --- последовательности, получаемые при секвенировании и содержащие информацию о фрагментах генома. Можно считать, что рид --- это строка четырех буквенного алфавита \textit{нуклеотидов}, соответствующая некоторому фрагменту стренда. При этом входные данные содержат ```вперемежку'' риды, соответствующие фрагментам различных стрендов (и обычно даже различных хромосом). %забиваем на colo space и N????

%\textit{Длина} рида --- количество нуклеотидов, его составляющих.

О риде, полученном по некоторому участку стренда говорят как о риде, \textit{покрывающем} данный участок. %участо кил фрагмент???

\textit{Покрытие} участка генома (в частности всего генома) --- среднее количество ридов, покрывающих отдельно взятый нуклеотид на этом участке. Обычно предполагают, что геном покрыт ридами более-менее равномерно.

Все секвенаторы совершают ошибки, наиболее частая из которых --- замена одного нуклеотида на другой (например, в риде находится T, а в соответствующей позиции генома --- A).
Для того чтобы помочь исследователям, в дальнейшем анализирующим эти данные, современные секвенаторы, вместе с каждым ридом, предоставляют \textit{характеристику качества} прочтения каждого из его нуклеотидов. Это информация о степени уверенности секвенатора в правильности того или иного нуклеотида.

Наиболее распространенной среди всех характеристик качества, являются так называемые Phred значения (Phred value coeffitient) $Q=-10log_{10}P$, где $P$ --- предполагаемая самим секвенатором вероятность ошибки в конкретном нуклеотиде.

К примеру, если Phred значение для некоторой позиции равно $30$, вероятность того, что на этом месте должен стоять другой нуклеотид, равна $1/1000$.

В силу технологических особенностей, качество нуклеотидов значительно снижается ближе к концу рида. И если для современных секвенаторов фирмы Illumina, Phred значения в начале и середине рида составляют, в среднем $30-40$, то к концу рида они падают до значений порядка $5$.

\textit{Расстоянием} между двумя участками одного стренда будем называть разность их начальных позиций.

Как уже отмечалось, большинство популярных секвенаторов имеют возможность генерации \textit{парных} ридов. Оба рида из пары соответствует участкам одного и того же стренда. Расстояние между этими участками будем называть \textit{внутренним расстоянием} парного рида.

Каждый набор парных ридов, соответствующих одному запуску секвенатора, характеризуется длиной одиночных ридов и предположительным внутренним расстоянием $d$.

В силу различных технологических причин, истинное внутреннее расстояние отдельно взятого парного рида может отличаться от $d$. 

Желаемое внутреннее расстояние $d$ производимых парных ридов задается при конкретном запуске секвенатора. О результатах такого запуска говорят как о множестве ридов (\textit{библиотеке ридов}) с предположительным внутренним расстоянием $d$. Ограничения на возможные значения $d$ варьируются от секвенатора к секвенатору и обычно находятся в диапозоне от $100$ до $10000$ нуклеотидов.

$k$\textit{-мер} --- строка длины $k$ в алфавите A,T,G,C. ($k$-mer$\in\{A,C,G,T\}^k$)

%\textit{Библиотека парных ридов с внутренним расстоянием $d$}  --- множество парных ридов

\subsection{Граф де Брюина}

Граф де Брюина  --- компактная структура, основанная на $k$-мерах, оказавшаяся крайне полезной для представления множества коротких ридов, обеспечивающих высокое покрытие графа.

Подход к сборке генома, основанный на использовании этих графов был предложен Павлом Певзнером с соавторами в 2000 году и с тех пор стал основным подходом к ассемблированию данных от NGS (Next Generation Sequencing) технологий. 

Более того, на сегодняшний день, это по сути единственный подход, которым можно хотя бы пытаться собирать большие геномы, содержащие миллиарды оснований, по коротким ридам $\approx 200$ нуклеотидов.

\subsubsection{Классическое определение}
В 1946 году, голландский математик Николас де Брюин занимался проблемой поиска минимальной циклической подстроки, содержащей все возможные двоичные последовательности длины $l$ в качестве подстрок.

Он рассматривал особый класс графов, определяемый следующим образом: рассмотрим $n-$буквенный алфавит и некоторое фиксированное число $l$. Сформируем множество всех возможных строк длины $l-1$ над этим алфавитом. Де Брюин рассмотрел граф $B(n, l)$, в котором в качестве вершин стоят все элементы построенного множества, а две вершины $v_1, v_2$ соединены направленным ребром тогда и только тогда, когда существует слово длины $l$, чьим префиксом является строка из $v_1$, а суффиксом строка из $v_2$. Граф содержит $n^{l-1}$ вершин. Нетрудно также заметить, что входящая и исходящая степень любой вершины этого графа равна $n$.

\subsubsection{Граф де Брюина для идеального множества ридов}
Построим подобную конструкцию по множеству ридов и выясним чему соответствует исходный геном в получившемся графе.

Но, для начала, упростим ситуацию, предположив, что множество наших ридов $S$ \textit{идеально}:
\begin{itemize}
\item Все риды имеют одну и ту же длину $l$.
\item Все возможные подстроки генома длины $l$ встречаются в множестве ридов.
\item Все риды соответствуют одному стренду (далее в этом разделе слова геном и стренд являются синонимами).
\item Риды были сгенерированы без ошибок.
\end{itemize}

Также предположим пока, что в восстанавливаемом геноме нет повторяющихся фрагментов длины $l$.

Граф де Брюина для такого множества ридов определим следующим образом.

Рассмотрим все различные префиксы и суффиксы ридов длины $l-1$. Построим граф $D$, в котором они являются вершинами и соединим две вершины $v_1$ и $v_2$ направленным ребром тогда и только тогда, когда существует рид с префиксом $v_1$ и суффиксом $v_2$. Каждый исходный рид соответствует ребру полученного графа.

Ясно, что получившийся граф $D$ является подграфом $B(4,l)$, так как получился тем же методом построения просто по меньшему множеству строк. 

Каждый путь в графе соответствует некоторой строке. Начинаем с $k-1$- мера, соответствующего первой вершине пути. Проходя по каждому ребру, приписываем последний нуклеотид соответствующего $k$-мера справа к имеющейся строке. В итоге получится строка длины $k-1+$ (количество ребер пути).

Нетрудно заметить, что в условии сделанных предположений, стренд соответствует некоторому \textit{эйлеровому} пути в графе $D$. Эйлеров путь --- путь, проходящий по каждому ребру ровно по одному разу. Также ясно, все эйлеровы пути в полученном соответствует строкам одинаковой длины, каждая из которых является решением задачи о кратчайшей надстроке (SSP, \cite{superstring}) для множества $S$.

\textbf{Замечание.} Для нахождения эйлерова пути в графе есть простой линейный алгоритм. Получается, что нам удалось решить задачу SSP за полиномиальное время! Но сделать это нам удалось лишь за счет очень сильного дополнительного условия на множество входных строк: оно содержит все подстроки длины $l$ исходной строки.

Попробуем для начала отбросить единственное предположение, которое мы сделали про сам геном: предположение об уникальности всех $l$-меров в геноме. Как уже отмечалось, в геноме присутствует большое количество повторяющихся участков, так что это условие заведомо не выполняется.

Если бы нам помимо множества ридов были бы известны кратности соответствующих $l$-меров в геноме, то мы бы просто добавить в граф ребра соответствующей кратности. Нетрудно заметить, что при этом исходный геном продолжал бы соответствовать некоторому эйлеровому пути в графе. На практике, получить  информацию о точной кратности конкретного $l$-мера в геноме оказывается невозможно. Таким образом, ввести в граф правильную кратность ребер не получится, поэтому вообще не будем использовать кратные ребра. Это, конечно, еще не повод отказываться от выбранного подхода, просто нужно понимать, что в полученном графе геном больше не соответствует эйлеровому пути.  

Если нет никакой дополнительной информации, то разумно вновь искать геном среди решений задачи SSP для множества $S$. Все они в точности соответствуют кратчайшим путям, проходящим по каждому из ребер графа не менее одного раза. Задача поиска таких путей называется \textit{Задачей Китайского Почтальона} (Chinese Postman Problem). Все ее решения являются кандидатами в геном, но у нас не хватает инфорамции, чтобы отдать предпочтение кому-либо из них. %Если у нас есть не одиночные, а парные риды, то мы можем пытаться отбросить из рассмотрения пути, не согласованные с внутренними расстояниями парных ридов. 

Есть один нюанс. В общем случае, задача Китайского почтальона тоже NP-полна! Если бы мы рассматривали парные риды и пытались отбросить из рассмотрения пути, не согласованные с их внутренними расстояниями, то это только усложнило бы задачу. Получилось, что даже максимально упростив ситуацию, мы вновь пришли к труднорешаемой задаче. 

Но на самом деле все не так плохо и полученный граф может принести нам пользу. Предположим, что мы нашли некоторый путь, который гарантированно присутствует в любом решении задачи китайского почтальона (например, таким является любой путь, внутренние вершины которого имеют исходящую степень $1$). Тогда строка, соответствующая этому пути гарантированно присутствует в искомом геноме.

\subsubsection{Граф де Брюина для множества ридов}

К сожалению, описанный выше способ построения графа по идеальному множеству ридов на практике не приведет ни к чему хорошему. Что естественно, ведь множество настоящих ридов обычно не удовлетворяет ни одному из условий ``идеальности'', введенных в прошлом разделе.

Начнем с того, что риды могут иметь имеют неодинаковую длину. Но даже если бы это было так, то невыполнение второго условия приведет к тому, что граф ``распадется'' на большое количество компонент связности. Действительно, отсутствие в множестве ридов некоторого $l$-мера $s$, присутствующего в геноме, приведет к ``разрыву'' в графе де Брюина пути, соответствующего любому фрагменту генома, содержащему $s$. 

Неожиданное решение этих проблем в том, чтобы извлечь из ридов, на чью малую длину мы постоянно сетовали, еще более короткие подстроки длины $k$, по которым уже построить граф де Брюина.

Итак, опишем построение финального варианта графа более формально.
Пусть дано множество ридов $S = {s_1, \dots , s_n}$. Через $S_k$ обозначим множество всех $k$-меров, встречающихся в качестве подстрок по крайней мере одного рида из множества $S$. Выберем некоторой $k$ и определим граф де Брюина $D(S_k)$ следующим образом: все элементы множества $S_{k-1}$ являются вершинами графа, $k-1$-мер $v$ соединен с $k-1$-мером $w$ направленным ребром, если $S_k$ содержит $k$-мер, у которого первые $k-1$ нуклеотид совпадают с $v$, а последние $k-1$ --- c $w$. Таким образом каждый $k$-мер из $S_k$ соответствует ребру графа $D(S_k)$.

Если строить граф таким образом, то связность нарушается только в том случае, если некоторый $k$-мер генома вовсе не встречается в ридах.
Выбрав достаточно маленькое $k$, можно снизить вероятность этого события практически до нуля. Но при этом не стоит забывать, что чем меньше $k$, тем больше путей, соответствующих различным участкам генома будут пересекаться в графе из-за наличия в них одинакового $k-1$-мера. Обычно выбирают $20<k<40$ нуклеотидов, конкретное значение прежде всего должно зависеть от длины ридов, а также кратности и равномерности покрытия.


Попробуем отказаться от предполоения о том, что риды принадлежат одному и тому же стренду. На практике, в большинстве случаев, в процессе получения ридов, секвенатор действительно ``перемешивает'' данные, полученные с различных стрендов. 
В силу сопряженной структуры ДНК, можно предварительно дополнить $S$ сопряженными (реверс-комплементарными) ридами.
Тогда построенный граф де Брюина будет представлять из себя как бы ``склейку'' графов, соответствующих каждому из двух стрендов. Теперь путь, соответствующий стренду не должен проходить по всем ребрам.

Опианный граф де Брюина обладает очень важным свойством, делающим его столь полезным для сборки геномов из коротких ридов, обеспечивающих высокое покрытие генома: его компактность. Действительно, количество вершин и ребер графа линейно относительно длины генома и не зависит от количества ридов. 

Теперь обсудим последнее из сделанных ранее предположений: безошибочность ридов.
На самом деле в ридах присутствует большое количество ошибок различных типов. Из-за них, в графе появляется большое количество вершин и ребер, не соответствующих никакому из участков генома. Их наличие, с одной стороны усложняет дальнейшую работу с графом, запутывая его структуру, а с другой --- может привести к появлению соответствующих ошибочных участков в итоговой сборке. Так как размер графа де Брюина зависит от количества ошибок в ридах, то оказывается, что он все-таки зависит от количества самих ридов. Но при использовании процедур исправления ошибок, влияние этой зависимости на размер графа не оказывается фатальным.

Теперь достаточно трудно формально описать свойства, которым должен удовлетворять путь, соответствующий некоторому стренду. 
И тем не менее, не все так плохо. Раньше ведь мы знали ограничения и все-равно не имели полиномиальных алгоритмов поиска путей им удовлетворяющим. 

Представим себе, что нам удалось удалить все ошибочные ребра из графа. Тогда, к примеру, как и раньше любой путь без ответвлений будет гарантированно соответствовать участку некоторого стренда.

\subsubsection{Сжатый граф де Брюина}
Альтернативным, но во многих смыслах более удобным способом представления той же информации, что и граф де Брюина, является \textit{сжатый граф де Брюина}. 

Сжатый граф де Брюина  --- это граф, получаемый из графа де Брюина в результате применения следующего преобразования сжатия неразветвленных участков:
\begin{algorithmic}
\WHILE{В графе есть вершина $v$, входящая и исходящая степени которой равны $1$}
\STATE Заменим ребра, инцидентные $v$ одним ребром, помеченным последовательностью, считанной по пути из удаляемых ребер %(конкатенация с учетом перекрытия на $k$)
\STATE Удалим $v$ из графа
\ENDWHILE 
\end{algorithmic}

Множество вершин полученного графа, является подмножеством вершин исходного графа ($k-1$- меров). Ребра соответствуют последовательностям нуклеотидов различной длины (не короче $k$). При этом последовательности на ребрах согласованы с последовательностями в вершинах (последовательность в ребре начинается и заканчивается с $k-1$-меров из соответствующих вершин). 

Последовательность, соответствующая пути в сжатом графе определяется также как и в несжатом, с той лишь разницей, что, проходя по ребру, мы добавляем в конец все нуклеотиды ребра, начиная с $k$-го.

Длиной ребра $e$ ($length(e)$) будем называть количество $k$-меров, содержащихся в этом ребре. То есть $length(e)=$ (количество нуклеотидов в ребре) $-(k-1)$. При таком определении, до сжатия графа все ребра имеют единичную длину.  

В дальнейшем, если не сделано специальных оговорок под графом де Брюина будет пониматься именно его сжатая версия.

\subsection{Общая схема работы ассемблера, основанного на графах де Брюина}

Общая схема работы современных ассемблеров, основанных на графах де Брюина:
\begin{itemize}
\item Предварительное исправление ошибок
\item Построение  графа де Брюина
\item Коррекция графа
\item Разрешение повторов
\item Получение контигов и скэффолдов
\item Этап консенсуса
\end{itemize}

Вкратце опишем что происходит на каждом из этих этапов.

\subsubsection{Предварительное исправление ошибок}

Как уже неоднократно отмечалось, в ридах присутствует большое количество ошибок. К примеру, для типичных данных, полученных секвенатором фирмы Illumina ошибки присутствуют в среднем в кажом третьем риде (если точнее, в 32\% ридов).
Если ничего не предпринять, то мы получим граф, содержащий большое количество ошибочных элементов. К примеру, в проекте по секвенированию \textit{N.meningitidis} (\cite{CP08}) несжатый граф, построенный по эмулированным данным состоял из 4,039,248 вершин, в то время как граф, построенный по реальным данным (для 20-меров), содержал 9,474,411 вершину. 

Все современные ассемблеры производят этап предварительной коррекции ошибок. 

Обычно, в начале, у ридов обрезаются концы низкого качества (на основе Phred коэффициентов).

Затем все подходы используют следующую идеею: каждая подстрока генома фиксированной длины (обычно, порядка 20-30 нуклеотидов), должна встречаться в ридах большое количество раз. Задача состоит в том, чтобы внести в риды в некотором смысле наименьшее (с учетом весовых PHRED коэффициентов) количество исправлений таким образом, чтобы все присутствующие в них подстроки, встречались бы в них достаточное (порог зависит от покрытия) количество раз.

Эти техники устраняют только ошибки вида ``замена одного нуклеотида на другой''. Но учитывая, для современных секвенаторов такие ошибки составляют абсолютное большинство, то процедуры устранения ошибок способны устранить до 98\% всех ошибок, находящихся в ридах. 

Так, в уже упомянутом проекте по секвенированию \textit{N.meningitidis} удалось снизить количество вершин до 4,081,857.

\subsubsection{Построение графа де Брюина}
Почти все популярные ассемблеры в явном виде строят некоторый сжатый граф де Брюина. 

\subsubsection{Коррекция графа}
Большая часть ошибок устраняется на этапе предварительного исправления (этап 1), но даже относительно небольшая часть ошибок, оставшаяся во входных данных приводит к тому, что количество ребер сжатого графа возрастает в десятки раз.

Для того, чтобы сделать результат сборки более качественным, предпринимаются попытки устранения из графа участков, соответствующих ошибкам входных данных.

%На этом этап этапе большую роль играет покрытие ребер, правильно

\subsubsection{Разрешение повторов}

Как уже было отмечено, в числе всех трудностей, значительно усложняющих задачу ассемблирования, особое значение играет наличие в геноме так называемых повторов. Это идентичные (или имеющие малое редакционное расстояние) участки генома, повторяющиеся 2 и более раз. Повторы имеют свою классификацию, которая не будет нас интересовать. Более 50\% длины генома млекопитающих занимают повторы различных типов. Это обусловлено биологическими причинами: не до конца изученными преобразованиями, за счет которых происходит эволюция генома.

В контексте графа де Брюина, построенного по $k$-мерам, \textit{повторами}, будем считать, все последовательности нуклеотидов длины $\ge k-1$, присутствующие в геноме более одного раза (быть может, на разных стрендах). Ясно, что каждый повтор в геноме ``запутывает'' графа де Брюина, склеивая участки графа, соответствующие различным участкам генома. Эти участки графа будем 

При разрешении повтора, происходит ``расклеивание'' некоторых ребер графа так, что из графа исчезают пути, которые не соответствуют никаким участкам генома.

Рассмотрим простой пример: представим себе два участка генома, содержащих повтор. Наличие повтора приводит к тому, что пути в сжатом графе, соответствующие этим участкам ``накладываются'' друг на друга (cм. рис.\ref{ris:repeat}а). Через ребро, соответствующее самому повтору при этом проходит уже 4 различных пути, лишь 2 из которых соответствуют исходным участкам. 

\begin{figure}[h]
\begin{minipage}[h]{0.49\linewidth}
\center{\includegraphics[width=0.5\linewidth]{repeat} \\ а)}
\end{minipage}
\hfill
\begin{minipage}[h]{0.49\linewidth}
\center{\includegraphics[width=0.5\linewidth]{repeat_resolved} \\ б)}
\end{minipage}
\caption{Иллюстрация разрешения повтора}
\label{ris:repeat}
\end{figure}

Разрешение повторов --- это этап, на котором происходит попытка обнаружения ``истинных'' путей в графе, проходящих через повторы, или, если немного переформулировать, связывание входов в повтор с выходами из него. Затем происходит разделение повтора, при котором создается по копии повтора для каждого истинного пути. (cм. рис.\ref{ris:repeat} б)

За счет чего можно отличить истинные пути? 

Вернемся к нашему примеру. Рассмотрим риды, проходящие через участок графа с рис. \ref{ris:repeat}а. Предположим, что есть риды, проходящие одновременно по обоим красным ребрам и есть риды, проходящие по обоим зеленым ребрам, в то время как ридов, идущих через красное и зеленое ребро не существует. Тогда мы можем предположить, что ``перекрестные'' пути являются побочными. 

Ясно, что подобные рассуждения могут принести плоды только в случае, если длина повтора меньше длины ридов (в противном случае по ним нельзя сделать вывод о связях между входами и выходами из повтора). Чем длиннее риды, тем больше повторов мы сможем разрешить. К сожалению, нам приходится работать с достаточно короткими ридами.

На помощь приходит то, что мы имеем дело с парными ридами. Вновь рассмотрим ситуацию на рис. \ref{ris:repeat}а, но в этот раз будем считать, что длина повтора превышает длину рида. При этом, если есть парные риды, ложащиеся обоими концами на красные ребра и есть --- ложащиеся на зеленые, но нет парных ридов, одним концом ложащихся на зеленое, а другим --- на красное ребро, то мы вновь можем разрешить этот повтор 

Общий итог такой: при использовании парных ридов у нас появляется возможность разрешить повторы, длина которых не превышает длину их внутреннего расстояния. %При этом в них, конечно, содержится меньше информации, чем в одиночных ридах длины равной внутреннему расстоянию.

Важно понимать, что в ходе разрешения повторов, граф перестает быть графом де Брюина, так как один и тот же $k$-мер теперь может присутствовать более графе более одного раза.

\subsubsection{Получение контигов и скэффолдов}

Контиги и скэффолды --- конечный результат сборки.

\textit{Контигами (contigs)} называют последовательности нуклеотидов, выдаваемые ассемблером в качестве участков генома.

Наиболее простая и популярная стратегия получения контигов по графу заключается в том, чтобы взять в качестве контигов последовательности, соответствующие отдельным ребрам графа.

При этом про структуру графа обычно полностью забывают и больше ее не используют.

\textit{Скэффолдами (scaffolds)} называют линейные структуры из объединения контигов.

Процесс получения скэффолдов включает в себя:
\begin{itemize}
\item установление кратности контигов
\item установление порядка \myignore{и взаимной ориентации} контигов 
\item вычисление длин зазоров между контигами
\end{itemize}

Установление кратности контигов происходит, в основном по информации о покрытии соответствующих участков графа.

Ключевую роль при получении скэффолдов играет информация о расстояниях, полученная из парных ридов.


\subsubsection{Этап консенсуса}
Это финальный этап сборки, на котором происходит прикладывание исходных ридов к полученным контигам и исправление ошибок, которые могли быть внесены на предыдущих этапах.